First-Order Transition to Incommensurate Phase with Broken Lattice Rotation 
Symmetry in Frustrated Heisenberg Model 
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We study a finite-temperature phase transition in the two-dimensional classical Heisenberg 
model on a triangular lattice with a ferromagnetic nearest-neighbor interaction Ji and an 
antiferromagnetic third-nearest-neighbor interaction J3 using a Monte Carlo method. Apart 
from a trivial degeneracy corresponding to 0(3) spin rotations, the ground state for J3 7^ has 
a threefold degeneracy corresponding to 120 degree lattice rotations. We find that this model 
exhibits a first-order phase transition with the breaking of the threefold symmetry when the 
interaction ratio is J3/J1 = —3. 
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The magnetic compound NiGa2S4 is a rare example 
of a two-dimensional triangular lattice antiferromagnet 
and stays in a spin-disordered state at low tempera- 
tures. Interesting behaviors of this material have been 
observed by Nakatsuji and coworkers; 1-4 ' the magnetic 
part of the specific heat exhibits an unusual double-peak 
structure and the magnetic susceptibility gradually in- 
creases without divergent behavior as the temperature 
decreases. They have estimated the Weiss temperature 
9w — —80 K from the susceptibility results, where this 
temperature corresponds to the high-temperature peak 
observed in the specific heat. Below the low-temperature 
peak, the specific heat shows a T 2 dependence and a 
short-range noncollinear order develops. Indeed, the de- 
velopment of incommensurate (IC) spin correlation has 
been confirmed by magnetic neutron scattering experi- 
ment. This observation was explained by the mean- field 
calculation of a Heisenberg model with a dominant anti- 
ferromagnetic (AF) third-nearest-neighbor interaction. 

The antiferromagnetic Heisenberg model on the tri- 
angular lattice has been studied exhaustively. Kawa- 
mura and Miyashita argued 5 ^ that a classical model with 
only an AF nearest-neighbor (NN) interaction exhibits 
a topological phase transition driven by the dissocia- 
tion of vortices at a finite temperature. One of the the- 
oretical scenarios 6 ^ proposed for NiGa 2 S4 is based on 
this mechanism, taking into account the effect of bi- 
quadratic interactions. Another scenario 7-9 ^ is the spin 
nematic ordering based on the S = 1 spin model with NN 
bilinear-biquadratic interactions. Although these theo- 
retical works predicted a number of interesting physical 
phenomena, they do not fully take into account the char- 
acteristic spatial structure of NiGa2S4, that causes the 
IC phase. The IC phase was observed in the model with 
the AF NN and AF sccond-NN interactions. 10) If the 
second- NN interaction is dominant, there are three IC 
structures of the ground state that can be transformed 
to each other by 120 degree lattice rotations. It is easy 
to see that a similar threefold degeneracy exists when 
the third-NN interaction is dominant. Because of this 
discrete symmetry, a finite temperature phase transition 
can take place without violating the Mermin and Wag- 



ner's theorem. 11 ^ Therefore, in order to understand finite 
temperature properties of NiGa2S4, it is necessary to dis- 
cuss the possibility of such a phase transition. 

In the present study, we investigate a two-dimensional 
classical Heisenberg model on a triangular lattice with 
a ferromagnetic NN interaction J\ (< 0) and an anti- 
ferromagnetic third-NN interaction J3 (> 0). The model 
Hamiltonian is given by 



W = Jr J2 



E 



(1) 
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where Si is the vector spin of unit length. The first sum 
runs over NN pairs of sites and the second sum runs over 
third-NN pairs. Here, we apply the periodic boundary 
condition and assume that J3 is dominant (|Js| > |Ji|). 
Indeed, the existence of the dominant interaction J3 
is supported not only by photoemission spectroscopy 12 ^ 
but also by first-principles calculation. 13 ) We neglect the 
second- NN interaction Ji , because the magnetic neutron 
scattering result suggests that | J2 1 is much smaller than 
|Ji| or I J3 1. 1 ) In addition, the above-mentioned three- 
fold degenerate structure of the ground state does not 
change even if we introduce a weak second-NN inter- 
action. Therefore, we can neglect J2 without changing 
essential physics. 

We first study the ground-state classical spin config- 
uration of the Hamiltonian cq. (1). If there is no NN 
interaction J 1; the ground-state spin configuration on 
each of the four sublattices, which has a double period of 
the lattice, is a 120 degree structure. Since J\ generates 
correlation between sublattices, the spin configuration of 
each sublattice is distorted from the 120 degree struc- 
ture. Hence, the ground state becomes the IC state. This 
IC phase is described by a spiral configuration with the 
wave vector k that minimizes the Fourier transform of 
the interactions. There are six such wave vectors k in- 
side the first Brillouin zone of the triangular lattice. The 
spiral configuration is given by 



Si = R cos k ■ Ti — I sin k ■ r , 



(2) 



where R and I are two arbitrary orthogonal unit vectors 
and Ti is the position of site i in the real space on the 
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Fig. 1. Schematic picture of spin configuration at k = ±(fc,0). 
Axis 1 is characterized by k, and axes 2 and 3 by k/2. There is 
a threefold degeneracy corresponding to 120 degree lattice rota- 
tions. The spin configuration on each of the four sublattices is 
close to, but not exactly, the 120 degree structure. 



triangular lattice. In the case of Heisenberg spins, the 
difference between k and — k can be absorbed in the def- 
initions of R and I. Thus, there are three distinct groups 
of states corresponding to k = ±(fc,0), ±(ifc, ^k), and 

±(ifc, —^k) in the IC phase. The schematic picture of 
the spin configuration at k = ±(fc,0) is shown in Fig. 
1. The spiral spin configuration along one of the three 
axes can be characterized by the wave number k, while 
others by k/2. These three state groups are well sepa- 
rated only for IC wave numbers. In other words, if the 
ordering is commensurate, i.e., if the ground state is a 
60 degree structure, the 120 degree spatial rotation can 
be achieved as a result of continuous spin rotation with 
the total energy being fixed. However, a technical diffi- 
culty characteristic of the IC nature of the ground state 
arises from the fact that the wave vector cannot take an 
arbitrary value in the reciprocal space when the system 
is finite. Thus, the wave vector k of the ground state de- 
pends on the system size in an irregular fashion. For the 
interaction ratio J3/ Ji , we use the value —3 in the follow- 
ing calculations. According to the original experiment, ^ 
■h/Ji = — 3 ~ —10. Therefore, the present value is the 
lowest, which is consistent with the experimental result. 
We choose the smallest value for a technical reason. For 
larger values, relatively large systems are required for 
stabilizing the structure that converges to the true IC 
ground state in the thermodynamic limit. A more sys- 
tematic study for various J3/J1 values will be published 
elsewhere. 14 - 1 In this case, the wave number of the IC spin 
configuration is k = \k\ = 1.92188 in the thermodynamic 
limit. These results indicate the possibility of the spon- 
taneous breaking of the threefold symmetry at a finite 
temperature. 

We perform classical Monte Carlo simulations based 
on the standard heat-bath method. Each run contains 
2 ~ 10 x 10 6 Monte Carlo steps per spin at each temper- 
ature. We make 40 ~ 64 independent runs for each size 
to evaluate the statistical errors. Throughout this letter, 
the Boltzmann constant is set to unity. 
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Fig. 2. Temperature dependences of the (a) specific heat C, (b) 
uniform susceptibility \, and (c) internal energy E. 



To discuss a finite-temperature phase transition, wc 
calculate the specific heat C, the uniform susceptibility 
X, and the internal energy E. C and x are defined by 

1 {E 2 } - (E) 2 



C = 



X = 



L 2 T 2 
1 (M 2 ) 



M 



(3) 



(4) 



(5) 



where L is the system size and (• • • ) indicates the ther- 
mal average. We show the results in Fig. 2. The spe- 
cific heat exhibits a single peak, which is narrower and 
higher for larger systems, indicating a phase transition. 
As mentioned above, the peak position depends on the 
system size irregularly owing to the compatibility of the 
ordering wave vector with the system size. Therefore, we 
cannot reliably estimate the transition temperature. We 
introduce a characteristic temperature T C (L) depending 
on the lattice size, which is defined as the peak position 
of the specific heat. The uniform susceptibility and in- 
ternal energy decrease rapidly at around T C (L). Clearly, 
the phase transition of the present model is different from 
that of the AF Heisenberg model on the triangular lattice 
with no third-NN interactions. 5 -* 

While Fig. 2 indicates the phase transition, we study 
the energy distribution to find the order of this transi- 
tion. If the system is at the first-order transition temper- 



3 



ST 




Fig. 3. Energy distributions P(E) near T C (L) for L = 60 (T/J 3 
0.4205), 72 {T/J 3 = 0.4262), and 84 (T/J 3 = 0.4240). 
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ature, the energy distribution P(E) should be bimodal. 
The distribution P(E) at T = T C (L) is shown in Fig. 3. 
Again, we cannot avoid the irregularity due to the in- 
commensurability and finitencss of the system. However, 
we can at least see a clear evidence of the bimodal energy 
distribution, i.e., the valley in the middle of the distri- 
bution deepens with increasing lattice size. Therefore, 
we consider that the phase transition is of the first or- 
der. The same results have been obtained irrespective of 
whether the initial spin configuration is set to that of the 
IC order or a random spin configuration. Although not 
shown in Fig. 3, the distribution P{E) is singly peaked 
for L = 36 and 48. This indicates that it takes con- 
siderably large systems to observe the first sign of the 
first-order transition. 

From the ground state properties of the present model, 
it is natural to expect the spontaneous breaking of the 
threefold symmetry. Based on this expectation, we cal- 
culate average bond energies along the three axes sepa- 
rately. To be specific, 

1 



£ ^ = Jj2 E s *' s i (/* = 1,2,3). (6) 



('J)nn II axis A* 

We sort the three averages in descending order and define 
Ei, Ei, and E 3 as 



Ei = (max{ei,£ 2 ,£3}), 

E 2 = (mid{£i,£ 2 ,£3}), and 

E 3 = (min{ei,£ 2 ,£3}}- 



(7) 



In Fig. 4(a), we present the temperature dependence of 
the average direction-specified bond energies for L = 60. 
Ei and Ei increase but E3 decreases below T C (L). This 
result implies that the threefold symmetry is broken and 
one of the three axes is selected. The existence of such a 
characteristic axis can also be seen in a snap shot of the 
spin configuration. To study the anomalous behavior of 
the average direction-specified bond energies more quan- 
titatively, we estimate the energy difference defined by 
AE = Ei — E3. The temperature dependence of AE 
is shown in Fig. 4(b). AE abruptly increases around 



Fig. 4. Temperature dependences of the (a) direction-specified 
bond energies Ei, E 2 , and E 3 for L = 60 and (b) AE (= E1-E3) 
for various sizes. 



T C (L) and the slope of AE seems to diverge as the lat- 
tice size increases. From these results, we conclude that 
the first-order phase transition is accompanied by the 
spontaneous breaking of the threefold symmetry. 

In the two-dimensional classical Heisenberg spin sys- 
tem, the true long-range order of spins is prohibited by 
the Mermin and Wagner's theorem. 11 ) Therefore, the oc- 
currence of the spontaneous symmetry breaking must be 
attributed to an order parameter other than bare spins. 
We consider the order parameter 



Q 



T> E QW, 



(8) 



Q(r) = si ■ s 2 + s 2 ■ s 3 - 2s 3 ■ si, 



(9) 

where Q(r) is defined on an upward elementary triangle. 
The corner sites 1, 2, and 3 are numbered counterclock- 
wise for a triangle. The sum of eq. (8) runs over all the 
upward elementary triangles. The quantity Q takes a fi- 
nite value in the ordered phase and in the disordered 
phase. The temperature dependence of {Q 2 } is shown in 
Fig. 5(a). In the disordered phase, the correlation length 
£q of Q(r) can be obtained from the structure factor 
SQ(k ) using the Ornstein-Zernike form 15 ) 



J_ S Q (0) 
\k \V S Q (k Q ) 



- 1, 



(10) 



S Q (ko) = ^E^M^> e ifc °- ( "\ (11) 



where is the ordering vector of Q(r) and k is any 
wave vector close to the ordering vector. Note that eq. 
(10) does not yield the correct values below T C (L). In Fig. 
5(b), we show the temperature dependence of the corre- 
lation length. The correlation length discontinuously in- 
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Fig. 5. Temperature dependences of the (a) (Q 2 ) and (b) corre- 
lation length £q/L. 



creases and becomes larger than the lattice size at around 
T C (L). Thus, the Q{r) degree of freedom freezes at T C (L). 

We consider implications of the present results on the 
experimental results of NiGa 2 S4. The important points 
in the experimental results are that the IC phase emerges 
at low temperatures and the specific heat has a broad 
double-peak structure. The present results are consis- 
tent with the incommensurability observed in the neu- 
tron scattering experiments. On the other hand, the spe- 
cific heat of the present model only exhibits a single di- 
vergent peak. It corresponds to the low-temperature peak 
(around 10 K) observed in the experiment, because the 
experimental estimate J3 = 30 K 16 -* is too small to be re- 
sponsible for the high-temperature peak at around 80 K. 
Considering the scenario based on the topological transi- 
tion, 5 ) we may observe two peaks, one corresponding to 
the first-order transition we discuss in the present arti- 
cle, and the other corresponding to the topological tran- 
sition. However, we have not observed any sign of the 
topological transition in the present case of J3/ Ji = —3; 
there is no clear sign for a first-order transition in the 
experiment, and, as mentioned above, the higher transi- 
tion temperature is too high. Therefore, it is not likely 
that all experimental evidences can be explained by the 
present model. The inclusion of various other terms, such 
as spin anisotropy 17 ) and magnetic field, 18 ) will be dis- 
cussed elsewhere. 14 ^ 

It may be suitable to make a few comments on the 
first-order phase transition in the present model. The 
three-state Potts model is the representative model with 
threefold symmetry breaking. 19 ) In two dimensions, this 
model has a second-order phase transition at a finite 
temperature, in a strong contrast to the present model, 
despite the same type of symmetry breaking. This is 
not very surprising considering that the order of the 
phase transition is determined by physics on short- 
length scales. In the present model, short-range physics 



is strongly affected by the IC nature of the magnetic or- 
dering. As we have discussed above, the IC ordering may 
be closely related to the first-order transition observed 
in the present study, because the commensurate struc- 
ture (i.e., the exactly 60 degree structure) cannot ener- 
getically separate the three degenerate thermodynamic 
states. 

To summarize, we have studied the two-dimensional 
classical Heisenbcrg model on the triangular lattice with 
the ferromagnetic nearest-neighbor interaction J\ and 
the antiferromagnetic third-nearest-ncighbor interaction 
J3. We have found that the finite temperature phase 
transition is of the first order, and that this transition 
is accompanied by the breaking of the threefold symme- 
try. The incommensurate nature of the ground state may 
be essential to the present type of phase transition. 
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